      program ubi
c
c     Given the centre of a sphere and its radius, this program will find 
c     the indices of particles that lie in this sphere in an evolved 
c     simulation. Then it finds where these particles are at an earlier
c     time (the input must specify the 2 simulation epochs required) and
c     calculates what cube is required to enclose them all at a time t=0. 
c     Ascii files of the initial and final particle positions are produced 
c     for use with ubi.sm.
c
c     The code also produce a file 'mask' which defines the high res
c     region more precisely.
c

c
      implicit none
      include 'mask.inc'
      integer npmax,ngmax
      parameter(npmax=82000000,ngmax=300000)
      integer np,i,j,k,count,part
      real x0(npmax),y0(npmax),z0(npmax)
      real xx0,yy0,zz0,rr0

      integer nq,nlev
      real rat
      integer iargc
      character*512 ubifile,mask_file,initpos,finalpos


      integer ix,iy,iz,isid
      integer ixc,iyc,izc
      real xc1,xc2,xc3,xc4
      real yc1,yc2,yc3,yc4
      real zc1,zc2,zc3,zc4
      

      real x(npmax),y(npmax),z(npmax)
      real a,lbox,mp,xcc,ycc,zcc,xcen,ycen,zcen,afin
      real xc,yc,zc,side,xmin,xmax,ymin,ymax,zmin,zmax
      real dxmax,dxmin,dymax,dymin,dzmax,dzmin,dx,dy,dz,rap,radius
      integer npgrid,icen
      real xi,yi,zi,rmin,adist,adist2
      logical inbox,insphere
      character*120 posfile
      real*8  sumx,sumy,sumz
      real fracmask
      real qxc,qyc,qzc,qsz
    
      external rap

      if (iargc().ne.1) stop 'usage: ubi.exe input_file '
      call getarg(1,ubifile)


      do i=1,imaskdim
        do j=1,imaskdim
          do k=1,imaskdim
               imask(k,j,i) = 0
          enddo
        enddo
      enddo



      open(22,file=ubifile,status='old')
      read(22,*) posfile
      read(22,*) mask_file
      read(22,*) initpos
      read(22,*) nq


      close(22)


      call read_rob_file(posfile,np,lbox,x0,y0,z0,xx0,yy0,zz0,rr0)
      npgrid = np**(0.333333333) + 1.e-5

      print*,'np     = ',np
      print*,'npgrid = ',npgrid
      print*,'lbox   = ',lbox

      

      print*,'---------------------------------------------------'




c     Now must find bounding box in the initial conditions (1/Z = 0).
c

c      call initpos(icen,npgrid,lbox,xi,yi,zi)
c      call initpos2(icen,npgrid,lbox,xi,yi,zi,x0,y0,z0)
      xc = x0(1)
      yc = y0(1)
      zc = z0(1)
      print*,'Original position of central particle ',xc,yc,zc



      dxmin = lbox
      dxmax = -lbox
      dymin = lbox
      dymax = -lbox
      dzmin = lbox
      dzmax = -lbox
      do i=1,np
         dx=rap(x0(i)-xc,lbox)
         dy=rap(y0(i)-yc,lbox)
         dz=rap(z0(i)-zc,lbox)
         dxmax = max(dxmax,dx)
         dxmin = min(dxmin,dx)
         dymax = max(dymax,dy)
         dymin = min(dymin,dy)
         dzmax = max(dzmax,dz)
         dzmin = min(dzmin,dz)
      enddo

 107  format (6f14.6,i8)
      print*,'dxmax,dxmin,dymax,dymin,dzmax,dzmin ',
     * dxmax,dxmin,dymax,dymin,dzmax,dzmin 

      xmax=rap(xc+dxmax-lbox/2,lbox)+lbox/2
      xmin=rap(xc+dxmin-lbox/2,lbox)+lbox/2
      ymax=rap(yc+dymax-lbox/2,lbox)+lbox/2
      ymin=rap(yc+dymin-lbox/2,lbox)+lbox/2
      zmax=rap(zc+dzmax-lbox/2,lbox)+lbox/2
      zmin=rap(zc+dzmin-lbox/2,lbox)+lbox/2

      print*,'xmax,xmin,ymax,ymin,zmax,zmin ',xmax,xmin,ymax,ymin,zmax,zmin 

      xc=rap(xc-lbox/2+(dxmax+dxmin)/2,lbox)+lbox/2
      yc=rap(yc-lbox/2+(dymax+dymin)/2,lbox)+lbox/2
      zc=rap(zc-lbox/2+(dzmax+dzmin)/2,lbox)+lbox/2
      side=max(dxmax-dxmin,dymax-dymin,dzmax-dzmin)

      print*,'Box selected: ',xc,yc,zc,' Size = ',side
      print*,'npgrid=',npgrid
 
      ixc = (2.*xc + side)/lbox * real(npgrid) + 1.e-5
      iyc = (2.*yc + side)/lbox * real(npgrid) + 1.e-5
      izc = (2.*zc + side)/lbox * real(npgrid) + 1.e-5

      print*,'Now growing box the mesh size '


      side = 1.25*side


      print*,'pen box selected ',xc,yc,zc,' Size = ',side

      rat = real(nq+2)/real(nq)
      nlev = int(log10(lbox/side)/log10(rat) + 0.5)
      print*,'Value of nlev = ',nlev

      side = lbox/rat**nlev





       
        print*,'Final box selected ',xc,yc,zc,side    


        print*,'In mask params ',xc/lbox,yc/lbox,zc/lbox,side/lbox

c----------------------------------------------------------------------
c       Uncomment for interactive use
c        
c        qxc=xc/lbox
c        qyc=yc/lbox
c        qzc=zc/lbox
c        qsz=side/lbox
c        print*,'Enter desired values for mask params'              !either (a)
c        read*,qxc,qyc,qzc,qsz
c----------------------------------------------------------------------
 
c  Uncomment if for example you want to run automatic script 
        print*,'Choosing parameters automatically '
        qxc = xc/lbox
        qyc = yc/lbox                                             ! or (b)
        qzc = zc/lbox
        qsz = side/lbox
c        qsz = 0.3262313
c-----------------------------------------------------------------------


        xc = qxc*lbox
        yc = qyc*lbox
        zc = qzc*lbox
        side=qsz*lbox

         
        print*,'Mask pars selected: ',xc/lbox,yc/lbox,zc/lbox,side/lbox

c    Now define the masked region within the box

      do i=1,np

         dx=rap(x0(i)-xc,lbox)
         dy=rap(y0(i)-yc,lbox)
         dz=rap(z0(i)-zc,lbox)

         ix = (dx/side+0.5)*imaskdim + 1
         iy = (dy/side+0.5)*imaskdim + 1
         iz = (dz/side+0.5)*imaskdim + 1 
         imask(ix,iy,iz) = -1
      enddo


      call create_mask(fracmask)
      

      open (31,file=mask_file,
     &                 form='unformatted',status='unknown')
      write (31) imaskdim,xc/lbox,yc/lbox,zc/lbox,side/lbox
      write (31) xx0,yy0,zz0,rr0
      write (31) imask
      close(31)

      write(6,*)' imaskdim = ',imaskdim

      stop
     




      print*,'Parameters of small box for pmdata_clus+pmdata_clus1 are:'
      print*,'xc,yc,zc,side ',xc/lbox,yc/lbox,zc/lbox,side/lbox
      print*,'bbox size Mpc/h ',side


      print*,'Mmask/Mblob = ', (576.0*side/lbox)**3*fracmask/real(count)



      close(18)

      end
c     
******************************************************************************
c
c     Function to do a periodic wrap-around in a box length l.
c
      function rap(diff,l)
      implicit none
      real rap,diff,l
      rap = diff
      if (diff .ge. l/2) then
         rap = rap - l
      endif
      if (diff .lt. -l/2) then
         rap = rap + l
      endif
      end
c
******************************************************************************
      logical function insphere(x,y,z,xc,yc,zc,lbox,radius,adist)
      implicit none
      real x,y,z,xc,yc,zc,lbox,radius
      real difx,dify,difz
      real adist,rap
      external rap
      difx=rap(x-xc,lbox)
      dify=rap(y-yc,lbox)
      difz=rap(z-zc,lbox)
      adist = difx*difx + dify*dify + difz*difz
      insphere=(adist.le.radius**2)
      return
      end
c
      logical function inbox(x,y,z,xc,yc,zc,size,lbox)
      implicit none
      real x,y,z,xc,yc,zc,size,lbox
      real dtiny
      real difx,dify,difz,rap
      external rap

      dtiny = lbox*1.e-6
      difx=abs(rap(x-xc+dtiny,lbox))
      dify=abs(rap(y-yc+dtiny,lbox))
      difz=abs(rap(z-zc+dtiny,lbox))
      inbox=(difx.le.size/2.).and.(dify.le.size/2.)
     * .and.(difz.le.size/2.)
      return
      end
c





      subroutine create_mask(fracmask)

      include 'mask.inc'


      cen = real(imaskdim+1)/2.
 
      print*,'Stage 1 of making the mask ...'     

      do icube=0,(imaskdim-1)/2

        ipaint = 100+icube
        imin = 1+icube
        imax = imaskdim-icube
        ir = (imax-imin)/2
        do i=imin,imax
         do j=imin,imax    
          do k=imin,imax
          irad = int(0.1+max(abs(real(i)-cen),abs(real(j)-cen),abs(real(k)-cen)))
             if (irad.eq.ir) then   ! boundary cube
                call test_cell(i,j,k,irad,ipaint,ival)
                if (ival.eq.1) imask(i,j,k)= ipaint
             endif
          enddo
         enddo
        enddo
      enddo



      print*,'Stage 2 ...'

      do i=1,imaskdim
       do j=1,imaskdim
        do k=1,imaskdim
          call test_cell2(k,j,i,ival)
          if (ival.eq.1) imask(k,j,i) = 10
        enddo
       enddo
      enddo


      print*,'Stage 3 ...'

      isum = 0
      do i=1,imaskdim
       do j=1,imaskdim
        do k=1,imaskdim
          if (imask(k,j,i).lt.11) then
               imask(k,j,i) = 1
               isum = isum+1
          else
               imask(k,j,i) = 0
          endif
        enddo
       enddo
      enddo

      fracmask = real(isum)/real(imaskdim)**3
      print*,'Fraction occupied ...',fracmask

      isel =0

 44   continue

      isel = isel+1
      print*,'isel = ',isel

      do i=1,imaskdim
       do j=1,imaskdim
        do k=1,imaskdim
          call test_cell3(k,j,i,isel,ival)
          if (ival.eq.1) imask(k,j,i) = isel+1
        enddo
       enddo
      enddo




      isum = 0
      do i=1,imaskdim
       do j=1,imaskdim
        do k=1,imaskdim
            if (imask(k,j,i).ne.0) then
                  isum = isum + 1
            endif
        enddo
       enddo
      enddo  


      fracmask = real(isum)/real(imaskdim)**3
      print*,'Fraction occupied ...',fracmask

      if (fracmask.lt.1.0) goto 44
      
 


      
       isum = 0
       do i=1,imaskdim
       do j=1,imaskdim
        do k=1,imaskdim
            if (imask(k,j,i).eq.1) isum=isum+1
        enddo
       enddo
      enddo       

      fracmask = real(isum)/real(imaskdim)**3
      

      return
      end



      subroutine test_cell(i,j,k,irad,ipaint,ival)
      include 'mask.inc'
      integer ival


      iqual1 = 0
      idqual2 = 1

      do ix=-2,2
       do iy=-2,2
        do iz=-2,2
           ixx = i+ix
           iyy = j+iy
           izz = k+iz

           if (((ixx.ge.1).and.(ixx.le.imaskdim)).and.
     &         ((iyy.ge.1).and.(iyy.le.imaskdim)).and.
     &         ((izz.ge.1).and.(izz.le.imaskdim))) then

                icell = imask(ixx,iyy,izz)
                if ((icell.eq.ipaint-1).or.(irad.eq.(imaskdim-1)/2)) iqual1 = 1
                if (icell.eq.-1) idqual2 = 0
           endif
        enddo
       enddo
      enddo

      if (imask(i,j,k).eq.-1) idqual2 = 0


     

      if ((iqual1.eq.1).and.(idqual2.eq.1)) then
               ival = 1
      else
               ival = 0
      endif




      return
      end
 
                
                            
      subroutine test_cell2(i,j,k,ival)
      include 'mask.inc'
      integer ival


      if (imask(i,j,k).lt.10) then
           ival = 0
           return
      endif


      icount = 0
      do ix=-3,3
       do iy=-3,3
        do iz=-3,3
           ixx = i+ix
           iyy = j+iy
           izz = k+iz

           if (((ixx.ge.1).and.(ixx.le.imaskdim)).and.
     &         ((iyy.ge.1).and.(iyy.le.imaskdim)).and.
     &         ((izz.ge.1).and.(izz.le.imaskdim))) then
                if (imask(ixx,iyy,izz).lt.10) then
                        icount=icount+1
                endif
           endif
        enddo
       enddo
      enddo


      if (icount.gt.100) then 
          ival = 1
      else
          ival = 0
      endif



      return
      end

      subroutine test_cell3(i,j,k,isel,ival)
      include 'mask.inc'
      integer ival


      if ((imask(i,j,k).gt.0).and.(imask(i,j,k).le.isel)) then
           ival = 0
           return
      endif


      icount = 0
      do ix=-2,2
       do iy=-2,2
        do iz=-2,2
           ixx = i+ix
           iyy = j+iy
           izz = k+iz

           if (((ixx.ge.1).and.(ixx.le.imaskdim)).and.
     &         ((iyy.ge.1).and.(iyy.le.imaskdim)).and.
     &         ((izz.ge.1).and.(izz.le.imaskdim))) then
                if (imask(ixx,iyy,izz).eq.isel) then
                        ival = 1
                        return
                endif
           endif
        enddo
       enddo
      enddo


      ival = 0




      return
      end


      subroutine read_rob_file(posfile,np,lbox4,x0,y0,z0,xx0,yy0,zz0,rr0)
      implicit none
      character*(*) posfile
      integer np
      real*4 lbox4
      real*8 lbox
      real*4 x0(*),y0(*),z0(*)
      real*4 xx0,yy0,zz0,rr0

      integer npart(6)
      real*8 zred
      integer i,j

      np = 0

      open(16,file=posfile,status='old',form='unformatted')
      read (16) npart
      print*,npart
      read (16) zred
      print*,'redshift =',zred
      read (16) lbox
      print*,'boxsize =',lbox

      ! rcrain -- I added this so the details of the Lagrangian region of the original 
      ! particle selection can be passed through, to eventually be picked up by GADGET
      ! and written into snapshots, so TRAPHIC can pick a suitable radial boundary at
      ! which to stop adding radiative sources.
      read (16) xx0
      print*,'OriginalSelection_x =',xx0
      read (16) yy0
      print*,'OriginalSelection_y =',yy0
      read (16) zz0
      print*,'OriginalSelection_z =',zz0
      read (16) rr0
      print*,'OriginalSelection_r =',rr0

      lbox4 = real(lbox)
      do j=1,6
       read (16) (x0(i+np),y0(i+np),z0(i+np),i=1,npart(j))
       np = np + npart(j)
      enddo
      close(16)


      end








